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Abstract 

This report presents preliminary results of an investigation into the development of a procedure 
to provide curvature continuity between biparametric cubic Bezier surface patches in the computer- 
aided design package known as SMART (Solid Modeling Aerospace Reasearch Tools). This initial 
effort was aimed at providing the designer with the ability to locally impose curvature continuity at the 
intersection of two Bezier curves without disrupting either the curvature or slope continuity that may 
exist at the ends of these curves. Such a method was found if the origianl Bezier control points are all 
coplanar. If they are not then it is possible to find a minimum deviation from exact curvature continuity. 
In cases where this is not sufficient, then an entire piecewise curve must be made curvature continuous 
simultaneously. A method was developed based on cubic splines which is very fast The procedure 
returns new Bezier control points which have both slope and curvature continuity. 

Introduction 

The development of computer-aided design packages have advanced to the state where it is now 
possible to build up a geometric description of an aerospace system in a very short time using high power 
graphics workstations. At the same time, it becomes desirable to have the ability to apply various 
numerical flow solvers to the resulting shapes to test their aerodynamic efficiency. These numerical 
flow solvers usually need the domain partitioned into small cells to provide a computational grid suitable 
for the discretized equations. 

The ultimate aim of this project is provide the ability of a computer-aided design package, 
SMART, developed at Langley Research Center, to include the generation of finite difference grids in 
a computational domain around the designed surfaces. Several methods for generating such grids are 
available but none are designed to smooth surface data. Hence if the surface has discontinuities the grid 
will be generated accordingly. The subsequent computation of the metrics of a coordinate transformation 
by using finite differences will thus have larger errors than usual. This inaccuracy of finite difference 
approximation to a differentioal equation will be larger than desired. To prevent these unwanted errors, 
it was necessary to begin with the requirement that the design surfaces must have curvature continuity 
over most of the domain. 
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Most computer-aided design packages describe the design surfaces by combining several surface 
elements called patches. These patches can be geometrically described by analytic functions of various 
kinds. Of particular interest are cubic polynomials. Here, the cubic polynomials are cast as biparametric 
Bezier curves. Cubic Bezier curves use control points as weights for a series of Bernstein polynomials and 
provide a relatively powerful mechanism for the designer to manipulate the shape, which is done by moving 
the control points. In the same manner, movement of the control points can allow for postprocessing of the 
designed shape as we wish to do here to impose curvature continuity. Thus, after the designercreates a shape, 
the resulting patches can then be smoothed by manipulating the location of the patch control points. 

Unfortunately, curvature continuity, along with slope continuity, requires information to be used 
from more than one patch. This couples all patches together requiring the simultaneous solution of new 
control point locations. If the number of patches is large, then this would require the solution of a rather large 
linear system and would likely be time-consuming even on high speed graphics workstations. It was thus 
desirable to investigate whether it would be possible to provide the curvature continuity for patches 
individually or at least with some minimal effort 

The examination of the problem took the form of first looking at how the designer actually constructs 
a shape and seeing if there were some way of providing a priori curvature continuity rather than at the end 
of the process. It seems that a significant amount of time is spent in the generation of two dimensional cross 
sections of the design shape. Thus cubic Bezier curves are generated describing these cross sections. 

Once the cross sections are described, they are assembled axially to build up the surface. A number 
of axial curves are then generated to connect the cross sections. These are also piecewise cubic curves with 
one cubic in between each cross section. In this manner, the surface is thus described by two families of 
piecewise cubic curves. The surface can be seen as a collection of four-sided patches with one cubic on each 
side. The task of providing curvature continuity between the patches is made much easier if there already 
exists curvature continuity in both families of curves. 

To obtain curvature continuity on these curves two tasks were identified: 

Task 1- Provide a means of adjusting the Bezier control points at the intersection of two cubic Bezier curves 
to gain slope and curvature continuity without disrupting the same at neighboring interserctions. 

Task 2 - Provide a means of ensuring curvature continuity on an entire series of cubic Bezier curves. 

Both of these tasks have been completed. The remainder of this report explains the approach, 
method, and presents an example for each task. The remaining effort of providing curvature continuity for 
entire network of patches will be the subject of the final report. 
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Task 1 - Curvature Continuity between adjacent Cubic Bezier Curves 
while maintaining Cl & C2 at the ends and Cl at the junction 

It is desirable, in CAD systems, when using cubic curves, to be able to attain various levels of 
continuity between adjacent curves. Here the desire is to attain C2 at intersections without disturbing 
the Cl and C2 continuity at neighboring intersections. In addition, it is assumed that the Cl continuity 
which already exists at the intersection will not be disturbed. By this is meant not only that the Bezier 
control points on either side of the intersction are colinear with it, but that the ratio of their distances from 
it will remain constant This ratio is also required at the neighboring junction points and does not allow 
for the control points nearest to neighboring junctions to be moved. It is thus allowed only to move, in 
space, the intersection point in question and the two control points adjacent to it 

If such a shift of control points is possible, this would allow the designer the flexibility of select- 
ing local junctures for enhanced continuity without having to simultaneously perform the operation at 
all junctures where such continuity is desired. In the present effort, it will be shown that if adjacent Bezier 
cubics have coplanar control points, there is a unique solution to the problem. If the control points are 
not all coplanar, then it is possible to find the minimum C2 discontinuity at the intersection. 

Analysis 

Consider the intersection of two cubic Bezier curves shown in Fig. 1. The control points are 
coordinate vectors denoted by b n , b n , , etc. Shown also are the requirements for C2 continuity at b n , 
ie. an auxiliary point d exists which not only is the intersection of extensions from the control polygon, 
but that ratios of certain lengths must also be held. 



Fig. 1 Adjacent Bezier curves, associated control points and 
auxilliary point d. The t and 1-t are meant to indicate ratios only. 
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When C2 continuity does not exist, there are two points, d, and dj, as shown in Fig. 2. This is the 
case more typically found. To gain C2 at b B , it is necessary to move one or more of the control 
points so that the two points d ( and are coincident 





Fig. 2 Adjacent Bezier curves when C2 at the juncture does not exist 
The auxilliary points dj and d^ are not coincident 


This can be done if the following are satisfied: 

b B . l = (l-t)b B _ 2 +td 
b B = (l-t)b n ., + tb n+1 
b B+t = (l-t)d + tb B+2 

To maintain C2 at the ends the following must also be satisfied: 

bo-3 * 2b B . 2 + b B .j = kj(b B .3 - 2b B . 2 + b B .j) 
b B+ j * 2b B+2 + b B+l = k 2 (b B+J - 2b B4 . 2 + b B+ , 

Where the overbar indicates the original position. Maintenance of Cl with distance ratios at the ends 
means that the points b o 2 and b n+2 can not be moved. Hence, only b B ,,b B , and b n+1 are to be relo- 
cated. Equations 1-3 further require that all but b n 3 and b B+3 must be coplanar for the possiblity of C2 
atb to exist 

The geometric interpretation of the problem is shown in Fig. 3. Equations 4 and 5 require 
that b^j and b B+1 can be moved only along lines parallel to the vectors represented by the second dif- 
ferences of the end control points. Thus kj controls the location of b^ and kj controls the loca- 
tion of b^ . Equations 1 and 3 then locate the auxilliary points dj and dj . Thus different values of 
k, move b^-and dj through parallel lines in space. Similarly, kj moves b^j and dj through two 
other parallel lines in space. It should now be obvious that if all control points are coplanar, then the 
lines swept by d, and dj are also coplanar and their intersection represents the solution of the prob- 
lem. Mathematically, we note that there are S equations for die 4 unknowns (the three control points 
and the auxilliar y point d) and two undetermined parameters k t and k 2 . This shows that if all of the 
control points for both curves are coplanar, then each point has two independent components. Hence 
the two parameters (kj and kj) are sufficient to make up a 5th unknown and the system is well-posed. 


( 1 ) 

( 2 ) 

(3) 


(4) 
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Fig. 3 Geometric interpretation of the solution of equations 1-5. As kj and kj 
are changed, the auxilliary points d and d move along lines parallel to the 
endpoint vectors. The intersection of these lines at d is the solution. 


Non-Coplanar Endpoints 

If b and b +3 are not coplanar with the other control points then the lines swept by dj and 
d. are skew"and no solution exists. However, there is an extremely useful interpretation of the 
closest approach of the two lines. The closest approach of the two skew lines represents the best 
possible match of the curvatures of the two cubic Bezier space curves. In fact, as the two endpoints 
approach being coplanar with the remaining control points, the two skew lines come closer to inter- 
secting and the respective curvatures of the two lines approach being exactly C2. This means that if 
the endpoints are not far from being coplanar with the rest of the control points then it is possible to 
attain near C2 continuity. 

To see why a unique solution may not exist when the endpoints are not coplanar, consider 
fitting two cubic splines between three points in space. Both Cl and C2 continuity are attained at 
the middle point but only specified slopes at the ends can be given. To attain specified curvatures at 
the ends, the middle point can be moved. There is a unique solution to this problem, even if the 
endpoint slopes are not coplanar. If they are, then this corresponds to the unique solution given in 
the last section. For the case when they are not, one need only examine the corresponding Bezier 
control points to see if those nearest the ends match the ones given in the problem. It is unlikely. 

We next return to the idea of finding this nearest C2 continuity. It amounts to finding the 
values of k. and kj for which the two skew lines have their closest approach. To find these values it 
is merely sufficient to find the location in space of the points on each line at closest approach. The 
following is an outline of the procedure and an example to show how near Cl continuity is attained. 





Procedure f or Nearest C2 
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Consider the skew lines in space lj and 1^ swept out by d t and d^ as the values of kj 
are changed as shown in Fig. 4. Note that kj and k 4 are directly proportional to kj and kj . 
equations for 1, and lj are given by: 

K =d i + K c x 

Ijsdj +k 4 c 2 


andk 4 

The 


• ( 6 ) 


Where c, and c 2 are the curvature vectors of the end points. Let the intersections of the line containing 
the shortest distance between the skew lines with the lines be called P and Q. We now seek to find P and 
Q and the associated values of k, and k 4 . 



Fig. 4 Skew lines given by sweep of dj and due to changing kj and k 4 . 


(7) 

L The cross product of the endpoint curvature vectors c, andc 2 defines the direction of 
the line containing P and Q. The distance between P and Q is found by projecting the line segment 
dj dj on the unit vector in the direction of PQ: 

w=(dA) te) <8) 

2. Next 1, is projected onto the plane perpendicular to PQ containing ^ . The equation 
for the projected line is: 

l^d, +k,c 1 + PQ (9) 

The projected line crosses lj atQ. Hence equating the projected line with lj gives an expression 
for solving for k, and k 4 . 

3. Having found kj and k 4 , b n I is found from equation 4, b n+1 is found from equation 
5, and b a is found from equation 2. Note that the original Cl and C2 at the ends and Cl at the inter- 
section are maintained exactly while C2 at the intersection is approximate. 
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Example (Coplanar contr ol points) 

To insure that a coplanar set of points was generated while trying to be general, a set of points 
was generated all with the same z-coordinate. The set was rotated about the x-axis and then about the 

y-axis. The resulting set of coordinates were used: 



n-3 

n-2 

n-1 

n 

n+1 

n+2 

n+3 

X 

y 

z 

-1.39545 

.76603 

1.58301 

-.71244 

1.63205 

1.76603 

.13878 

2.48076 

1.84038 

.88660 

2.67128 

1.53564 

2.00833 

2.95707 

1.07853 

■ 

3.08109 

.07321 

-1.46340 


The mismatch in initial curvature at the intersection is shown best by noting that equations 1- 
3 can also be expressed as requiring: 

b. - 2b.., + b B _ 2 = k(b„ - 2b n+l + b B+2 ) 

Hence each component of the second difference vector on one side of the intersection must be the same 
multiple (fraction) of the components on the other side. This means that the ratios of the components 

must be the same. That is: _ 2b n _i + b n _2 _ . ( 11 ) 

bn- 2b n+l + b n+2 


In this example the initial coordinates gives: 



k 

X 

y 

z 

3.2324 

.4086 

.3592 


Table I. Second difference ratios of coordinates at the intersection 
Thus the ratio of the second differences of the coordinates on either side of the junction are rather 



k 

X 

y 

z 

,4444 

AAAA 

• T V V V 

,4444 


Table n. Second difference ratios ot tne coordinates at the intersection 
after imposing curvature continuity 


Hence for the case of coplanar control points, the C2 is attained exactly. 


Example (No n coplanar) 

In this example, the first point, b n _ 3 , was displaced by a distance of .8 in the z-direction before 
the x- and y-axis rotations. This guaranteed that all points were coplanar except the first which had co- 
ordinates of (-1.04904, .36603, 2.18301)and was 35.6 deg out of the plane containing the other control 
points. The initial curvatures were the same as the above example, since the end control point is not 
included in computing the curvatures. The final curvatures, however, are given in Table HI below. 
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k 

X 

y 

z 

.3963 

.3996 

.3904 


Table HI. Second Difference ratios for noncopanra control points 
after curvature continuity is imposed. 

points. 

Hence the procedure outlined above allows for the designer to go to a single ™*™^™°** 

composite cubic Bezier curve and do cOTt^acoplMwrcoMrol points, this 

points of the two ^«"e co pan . S inttrsecdo ns where the control points on either side are 

ZZ+X Sows to come as close as possible. A program listing is given tn 

Appendix A which carries out the steps given above. 




Task 2- Curvature continuity (C2) of a series of cubic Bezier 
curves which are at first only CO 


9 


As it has been shown, if the initial control points are not coplanar, then exact curvature continuity, 
is locally unobtainable, though it is possible to generate the coefficients locally of that pair for Bezier curves 
which are as close as possible to curvature continuous. This may suffice for many applications. But if not, 
it be comes necessary to attain the continuity for all cubics in the piecewise curve simultaneously. 

Fortunately, this can be done relatively easily. The problem is very similar to setting up an 
interpolant between points in space. The procedure is one where the endpoints of each Bezier curve is kept 
fixed in space. Thus the outer Bezier control points of each cubic are retained. We will thus show that the 
problem is solved by replacing the interior control points with those computed from the corresponding cubic 
splines placed through the end points of each piecewise cubic. The end conditions used are thosewhich keep 
the original slopes at the ends of the piecewise curve. 


Analysis 


Consider the series of Beziercontrol points, control polygons, and corresponding cubic curves shown 
in Fig. 5 below. . b^ b JS 

ft © 


b„ b. 


21 a 

b„ 


boT- 


bo, ..*> 

/ b,r s k © * I, 

•' b 22 ->. b(n / b 34 

501 N b * 3 

® ® \ b B 

b M © 


y b. 


Rg. 5. Example of piecewise cubic Bezier curve where neither Cl nor C2 exist. 


For clarity, the control polygon has not been drawn. Since the Bezier control points on either side of the three 
interior intersections are not collinear with the control point at that intersection, it is clear that slope contiuity 
is absent Table IV shows the ratio of curvature vectors for each of the three interior intersections. As they 
are not equal, it is clear that curvature continuity is also absent. 


mm 

Intersection J 

■ 




D 

1.0769 

.8000 1.5000 

-2.6364 

H 

-6.4000 

-.1250 .7273 

.1250 


Table IV. Curvature ratios at the interior points of initial 
set of cubic Bezier curves. 

Now e nnyiriw the cubic spline which has been placed through the inersection points as shown in Rg. 6. The 
intersection points are also the outer two control points of any one cubic. Hence, only new inner control 
points need be found. To place a cubic spline through these points any standard cubic spline routine can 
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Fig. 6 Cubic spline through the intersections of the Bezier curves. 

be used. The user only need specify two end conditions. These end conditions can either be those called 
" natur al" wherein the second derivative is zero; "slope", where the slope at the end is specified; or some 
combination. It seemed most appropriate in this circumstance to use the slopes existing at the ends of original 
Bezier cubics. With these end conditions, a piecewise parametric cubic spline was passed through the ends 
and the in terse tion points. 

For this effort, an optimized parametric cubic spline routine was written where a single parameter, 
t, varied between zero and one on each cubic. Each cubic has the form: (12) 

x = a i t 3 + (5 i t 2 + Yit + *i 


where x represents the coordinates of a point on the cubic, and x i contains the coordinates of the ith cubic 
when t=0. Thus, if there are N cubics, there will be 3N unknowns. To get 3N equations we note that the 
condition of CO continuity requires; 



ai-i+Pi-i + Yi-i + Xi-^Xj 

i e [1,N] 

(13) 

Cl requires; 

3a w + 2p M +y i .,*Yi 

ie[l,N-l] 

(14) 

and C2 requires: 

30,-, + Pm -ft 

i€[l,N-l] 

(15) 


There are 3N equations and this leaves the two end conditions to get the remaining two equations. The usual 
practice is to use these equations to write a relationship between one of the coefficients at three different 
points. This eliminates the other two variables. Here we used eq. (13M15) to write and eqation for p: 



Pi-i + *Pi + Pi*i = 3( x «-i ~ 2x i + x i+0 i € [2,N-1] 

(16) 

The linear system is computed by adding the end conditions in terms of P . Hence we have: 



■^Pl + jP2~ ^{ X 2 X l) S 1 

(17) 


■jPn-i + = 3(x n+1 - x N ) - 2(x n - x N _i) - s 2 

(18) 

where 

Sj = 3(b u - b 01 ) 

(19) 


= 3(b 3N — b lN ) 



Eq uati ons (16)-(19) form a tridiagonal system of equations for each coordinate direction and can thus be 
solved very quickly. With each of the p's found, the other coefficients are computed from eq. (14) and (15). 
This defines the cubics in the piecewise curve. 
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Once the cubic splines have been generated, the inner two Bezier control points for each cubic can 
be found from: (20) 

•>2i = l>3i-|(3“i + 2 Pi + ri) 

The resulting new curves and corresponding control points are shown in Fig. 7. The dashed lines are 
the original cubics and the open circles are the original Bezier control points. It may be noted, as a visual 
check of slope continuity, that at each intersection, the two control points on either side of the intersection 

© 


Fig. 7 New cubic Bezier curve with new control points. 

point are now collinear, unlike the original control points. The control points next to those at the very ends 
of the total curve are the same as the original ones as a consequence of retaining the original slope at the end 
points. The fit of the new curve is close to the old except where the old curve was rather badly matched in 
slope and in curvature (ie. at intersection ). 

The new curvature ratios, shown below inTable V., show that curvature continuity has been 
achieved. 


Table V. Curvature ratios at the 4 interior intersections 

This completes the second task. The procedure allows for virtually any combination of cubic Bezier 
curves which describe a piecewise curve to be made into one which has both slope and curvature continuity. 
The original control points are retained at each of the curve intersections and new ones are created for the 
interior of each curve. The new control points are computed from the coefficients of the associated spline 
curve through theoriginal intersections. The control points adjacent to the outside ends of the first and last 
cubic have also been retained since the end requirement of slope matching the original curve has been 
im po sed It may be noted that there is no requirement on the original control points to be coplan ar as in 
task one. A program listing is given in Appendix B for a subroutine which takes the control points of a series 
of Bezier curves, fits a cubic spline through the intersections and endpoints, and recalculates the interior 

Bezier control points. 
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Appendix A 

Program Listing for Task 1 

The following program in FORTRAN will compute the values of kj and k* given an initial 
set of control points, and adjust the position of the control points b o ., , b .and b„ +I . 

SUBROUTINE C2(X,Y,Z) 

DIMENSION X(4, 2),Y(4, 2), Z(4, 2) 

Q. find curvatures of ends 

CEX1 = X(3,l) - 2*X(2,1) + X(l,l) 

CEY1 = Y(3,l) - 2*Y(2,1) + Y(l,l) 

CEZ1 = Z(3,l) - 2*Z(2,1) + Z0.1) 

CEX2 = X(4,2) - 2*X(3,2) + X(2,2) 

CEY2 = Y(4,2) - 2*Y(3,2) + Y(2,2) 

CFZ2 = Z(4,2) - 2*Z(3,2) + Z(2,2) 

C . FIND t TO MAINTAIN C2 AT b B 

DX = X(4,l) - X(3,l) 

DY = Y(4,l) - Y(3,l) 

DZ = Z(4,l) - Z(3,l) 

51 = SQRT(DX**2 + DT**2 + DZ**2) 

DX = X(2,2) - X(3,l) 

DY = Y(2,2) - Y(3,l) 

DZ = Z(2,2) - Z(3,l) 

52 = SQRT(DX**2 + DY**2 + DZ**2) 

T = S1/S2 


q FIND dj and dj 

XD1 = X(2,l) + (X(3,l) - X(2,l)) / T 
YD1 = Y(2,l) + (Y(3,l) - Y(2,l)) / T 
ZD1 = Z(2,l) + (Z(3,l) - Z(2,l» / T 

XD2 = (X(22) - T * X(3,2)) / (1 - T) 

YD2 = (Y(2 a) - T * Y(3,2)) / (1 - T) 

ZD2 = (Z(2^) - T * Z(3J2)) / (1 - T) 


C 


FIND P AND Q 


DX = XD2 - XD1 



DY = YD2 - YD1 
DZ = ZD2 - ZD1 
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CEX = CEY1 * CEZ2 - CEY2 * CEZ1 
CEY = CEX2 * CEZ1 - CEX1 * CEZ2 
CEZ = CEX1 * CEY2 - CEX2 * CEY1 

CL = SQRT ( CEX* *2 + CEY**2 + CEZ**2) 

PQ = ( DX * CEX + D Y * DE Y + DZ * CEZ) / CL 

CEX = CEX * PQ/CL 
CEY = CEY * PQ / CL 
CEZ = CEZ* PQ/CL 

DEN = CEX1 * CEY2 - CEX2 * CEY1 

K3 = ((CEY - DY) * CEX2 - (CEX - DX) * CEY2) /DEN 

K4 = ((CEY - DY) * CEX1 - (CEX - DX) * CEY1) /DEN 

PX = XD1 + K3 * CEX1 
PY = YD1 + K3 * CEY1 
PZ = ZD1 + K3 * CEZ1 

QX = XD2 + K4 * CEX2 
QY = YD2 + K4 * CEY2 
QZ = ZD2 + K4 * CEZ2 

NEW K1 AND K2 — 

K1 = (T * PX - (1 + T) * X(2,l) + X(l,l)) / CEX1 
K2 = ((1 - T) * QX + (T - 2) * X(3,2) + X(4,2)) / CEX2 

FIND NEW b , b , andb 

X(3,l) = 2 * X(2,l) • X(l,l) + K1 * CEX1 
Y(3,l) = 2 * Y(2,l) - Y(l,l) + K1 * CEY1 
Z(3,l) = 2 * Z(2,l) - Z(l t l) + K1 * CEZ1 

X(2J2) = 2 * X(3,2) - X(4,2) + K2 * CEX2 
Y(2^) = 2 * Y(3,2) - Y (4,2) + K2 * CEY2 
Z(2,2) = 2 * ZQ2) - Z(4,2) + K2 * CEZ2 

XN = (1 - T) * X(3,l) + T * X(2^) 

YN = (1 - T) * Y(3,l) + T * Y(2^) 

ZN = (1 - T) * Z(3,l) + T * Z(2,2) 

X(4,l) = XN 
Y(4,l) = YN 



Z(4,l) = ZN 


X(U) = XN 
Y(U) = YN 
Z(U) = ZN 


C 


THAT'S ALL 


RETURN 

END 
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Appendix B. 

Curvature Continuity on a Piecewise Cubic Curve 


The following FORTRAN subroutine is designed to perform the steps listed in Task 2. That is, given 

an input set of N Bezier cubic control points, cubic splines are put through the intersection points and new 
interior control points are found. 

SUBROUTINE BCRV(X,Y,Z,N) 


DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 


X(4,N),Y(4,N),Z(4,N) 

XX(N+1),YY(N+1)ZZ(N+1) 

A(N),B(N),C(N)D(N) 

AA(N),BB(N),CC(N) 


C- 

C 

c 

c 

c 

c 

c 

c 

c 


SOME OF THE VARIABLES 

N Number of Bezier Curves 

X,Y,Z Coordinates of Bezier Control points 

XX.YY.ZZ Coordinates of the ends and intersections of the piecewise curve 
A,B,C,D Coefficients of the cubic spline matrix 
AA.BB.CC Coefficients of the cubic splines 

FIND SLOPES AT THE ENDS 


SIX = 3. * (X(2,l) - X(l,l)) 
S1Y = 3. * (Y(2,l) - Y(l,l)) 
S1Z = 3. * (Z(2,l) - Z(l,l)) 


S2X = 3. * (X(4,N) - X(3,N» 
S2Y = 3. * (Y(4,N) - Y(3,N)) 
S2Z = 3. * (Z(4,N) - Z(3,N)) 


GET THE INTERSECTION AND END POINTS 

DO 1001= 1,N 


XX(I) = X(1,I) 
YY(I) = Y(1J) 
ZZ(I)=Z(U) 


100 CONTINUE 

XX(N+1) = X(4,N) 
YY(N+1) = Y(4N) 
ZZ(N+1) = Z(4,N) 



c . SET UP CUBIC SPLINE MATRIX 

C- LOOP ON DIMENSIONS 

DO 1000 K - 1,3 


DO 200 1 = 2.N-1 


C(D= 1 

B(I) = 4 
A(I) = 1 

IF (ICEQ.l) THEN 

dd = xxa+D - 2. * xx(d + xxa-D 

ELSEIF (K.EQ.2) THEN 
DD = YY(I+1) - 2. * YY(I) + YY(I-1) 
ELSH 

dd = zza+D - 2. * zza) + zza-D 

ENDIF 

D(I) = 3. * DD 
200 CONTINUE 


B(l) = 2. / 3. 

A(l)= 1 ./3. 

IF(K.EQ.1)THEN 
D(l) = 3. * (XX(2) - XX(1» - SIX 
ELSEIF (K.EQ.2) THEN 
D(l) = 3. * (YY(2) - YY(1)) - S1Y 
ELSE 

D(l) = 3. * (ZZ(2) - ZZ(1)) - S1Z 
ENDIF 


B(N) = 7./ 3. 

C(N) = 2. / 3. 

IF (ICEQ. 1) THEN 

D(N) = 3. * (XX(N+1) - XX (N)) - 2. * (XX(N) - XX(N-l)) - S2X 
ELSEIF (K.EQ.2) THEN 

D(N) = 3. * (YY(N+1) - YY(N)) - 2. * (YY(N) - YY(N-1)) - S2Y 
ELSE 

D(N) *= 3. * (ZZ(N+1) - ZZ(N)) - 2. * (ZZ(N) - ZZ(N-1)) - S2Z 
ENDIF 

C- SOLVE THE MATRIX 


CALL TSOLV(A,B,CJ),BB) 



c- — 


FIND THE OTHER CUBIC SPLINE COEFFICIENTS 


IF (K.EQ.1) THEN 


DO 300 1 = l.N-l 

AA(I) = (BB(I+l)-B(I))/3. 

CC(I) = xxa+l) - XX(I) - AA(I) - BB(I) 

300 CONTINUE 

CC(N) = 3. * AA(N-1) + 2. * BB(N-l) + CC(N-1) 
AA(N) = XX(N+1) - XX (N) - BB(N) - CC(N) 

ELSEIF (K.EQ.2) THEN 

DO 320 1 = l.N-1 
AA(I) = (BB(I+1) * BB(I)) / 3. 

CC(I) = YY(I+1) - YY(I) - AA(I) - BB(I) 

320 CONTINUE 

CC(N) = 3 *AA(N-1) + 2 *BB(N-1) + CC(N-1) 
AA(N) = YY(N+1) - YY(N) - BB(N) - CC(N) 

ELSE 


DO 3401 = l.N-l 

AA(I) = ( BB(I+1) - BB(I» / 3. 

CC(I) = ZZa+1) - ZZ(I) - AA(I) - BB(I) 

340 CONTINUE 

CC(N) = 3 *AA(N-1) + 2 *BB(N-1) + CC(N-l) 
AA(N) = ZZ(N+1) - ZZ(N) - BB(N) - CC(N) 

ENDIF 


NOW GET NEW BEZIER COEFFICIENTS - 

DO 500 1 = 1 TO N 
IF (KEQ.l) THEN 
DX0 = CC(I) 

DX1 =3. * AA(I) + 2. * BB(I) + CC(I) 

X(2,I) = X(U) + DX0/3. 

X(3J) = X(4J)-DXl/3. 

ELSEIF (K.EQ.2) THEN 

DY0 = CC(I) 

DY1 = 3. * AA(I) + 2. * BB(I) + CC(I) 



Y(2,I) = Y(U) + DY0/3. 

Y(34) = Y(4J) - DY1 / 3. 

ELSE 

DZO = CC(I) 

DZ1 = 3. * AA(I) + 2. * BB(I) + CC(I) 
2X2,1) = Z(1J) + DZO / 3. 

Z(3,I) = ZX(4J) - DZ1 / 3. 

ENDIF 

500 CONTINUE 

c END OF DIMENSION LOOP 

1000 CONTINUE 

C- THATS ALL 


STOP 

END 

C 

SUBROUTINE TSOLV(A,B,C,D,N,BB) 

DIMENSION A(N),B(N),C(N)T>(N)3B(N) 

C ELIMINATE C’S 

DO 1001 = 2, N 

CBI = C(I) /B(I-l) 

B(I) = B(I) - CBI * A(I-1) 
Da)=D(i)-cBi*Da-i) 

100 CONTINUE 

BB(N) = D(N) / B(N) 

DO200I = N-l,l,-l 

BB(I)' = (D(I) - A(I)*BB(I+1)) / B(I) 

200 CONTINUE 


RETURN 

END 




